CoCalc Logo Icon
DocsShareSupport News Sign UpSign In
Views: 100
Image: ubuntu2004
Embed | Download | Raw |
Kernel: Python 3 (ipykernel)
from __future__ import annotations import typing import random from IPython.core.display import SVG import pyomo.environ as pyo import pyomo.util.infeasible from pysat.solvers import Solver from pysat.formula import CNF import py_svg_combinatorics as psc from ipywidgets import widgets, HBox from collections import Counter from pprint import pprint from random import randint import numpy as np from IPython.display import IFrame import IPython from copy import copy import os from pathlib import Path nbname = '' try: nbname = __vsc_ipynb_file__ except: if 'COCALC_JUPYTER_FILENAME' in os.environ: nbname = os.environ['COCALC_JUPYTER_FILENAME'] title_ = Path(nbname).stem.replace('-', '_').title() IFrame(f'https://discopal.ispras.ru/index.php?title=Hardprob/{title_}&useskin=cleanmonobook', width=1280, height=300)

Генератор тестовых данных

def gen_random_inputs(k: int, n_vars: int, n_terms: int) -> np.ndarray['n_terms,k']: """ Returns C - the matrix of the CNF ter A var is an integer in [1, n_vars] A literal is a var, -var or 0, which means False for disjunctions (i.e. no effect) The set of all vars is implicitly defined as np.arange(n_vars) + 1 """ C: np.ndarray['n_terms,k'] = np.random.randint(n_vars + 1, size=(n_terms, k)) C *= np.random.choice([-1, 1], size=(n_terms, 3)) return C
k = 3 n_vars = 10 n_terms = 50 C = gen_random_inputs(k, n_vars, n_terms) k, n_vars, n_terms, C
(3, 10, 50, array([[ 0, -4, 6], [ 0, -8, -1], [ -8, 8, -2], [ 0, -3, -9], [ 1, -4, -8], [ 10, -3, 5], [-10, -8, 3], [ 7, -3, -8], [ -1, -10, -1], [ 1, -1, 8], [ -8, -5, 0], [ -5, 3, -4], [ -4, -10, 1], [ -2, 2, 2], [ 3, 0, 5], [ 5, -4, -6], [ -5, 4, 6], [ 10, 4, -9], [ -4, 8, -1], [ 7, -1, -9], [ -1, -9, -2], [ -5, -3, 0], [ 6, -9, -9], [ -2, 2, -8], [ -7, -3, -3], [ 8, 1, 0], [ 2, 9, 0], [ -4, -8, 2], [ 2, -1, 9], [ -2, 3, 5], [ -3, -9, 6], [ -9, 9, 0], [ -7, 1, 5], [ 4, 0, 0], [ 1, -9, -1], [ 9, -8, 1], [ 2, -2, 4], [ -3, -7, -10], [ -7, -10, -1], [ 10, 7, 6], [ 3, 4, 0], [ 7, -7, 4], [ 6, -2, -2], [ 3, 0, 8], [ -4, 5, -5], [ 10, 3, 0], [ -1, -7, -2], [ 3, 9, 4], [ 3, 4, -7], [ -8, 0, 8]]))

Pyomo-модель

def get_model(k: int, n_vars: int, n_terms: int, C: np.ndarray['n_terms,k']) -> pyo.ConcreteModel: m = pyo.ConcreteModel(name="maximum k-satisfiability") m.n_vars = n_vars m.n_terms = n_terms m.k = k """ C holds the matrix of the CNF terms. Each row is a term, each column is a literal. A literal is a var, -var or 0, which means False for disjunctions (i.e. no effect) A var is an integer in [1, n_vars] """ assert C.shape == (n_terms, k) m.C = C """ var_values holds the truth values of the variables. """ m.I = range(m.n_vars) m.J = range(m.n_terms) m.K = range(m.k) m.var_values = pyo.Var(m.I, domain=pyo.Binary) """ term_items holds the values of the literals in each term. """ m.term_items = pyo.Var(m.J, m.K, domain=pyo.Binary) @m.Constraint( m.J, m.K ) def evaluate_term_items(m: pyo.ConcreteModel, j: int, k: int) -> bool: if m.C[j, k] == 0: return m.term_items[j, k] == 0 if m.C[j, k] > 0: return m.term_items[j, k] == m.var_values[m.C[j, k] - 1] if m.C[j, k] < 0: return m.term_items[j, k] == 1 - m.var_values[-m.C[j, k] - 1] assert False, "Unreachable" """ term_values holds the values of the terms. Note: every term is a disjunction of literals """ m.term_values = pyo.Var(m.J, domain=pyo.Binary) @m.Constraint(m.J, m.K) def evaluate_terms_low(m: pyo.ConcreteModel, j: int, k: int) -> bool: return m.term_values[j] >= m.term_items[j, k] @m.Constraint(m.J) def evaluate_terms_high(m: pyo.ConcreteModel, j: int) -> bool: return m.term_values[j] <= sum(m.term_items[j, ...]) """ The objective is to maximize the number of satisfied terms. """ m.obj = pyo.Objective(expr=sum(m.term_values[...]), sense=pyo.minimize) return m m: pyo.ConcreteModel = get_model(k, n_vars, n_terms, C)
solver = pyo.SolverFactory('cbc') solver.solve(m).write()
# ========================================================== # = Solver Results = # ========================================================== # ---------------------------------------------------------- # Problem Information # ---------------------------------------------------------- Problem: - Name: unknown Lower bound: 35.0 Upper bound: 35.0 Number of objectives: 1 Number of constraints: 170 Number of variables: 60 Number of binary variables: 210 Number of integer variables: 210 Number of nonzeros: 50 Sense: minimize # ---------------------------------------------------------- # Solver Information # ---------------------------------------------------------- Solver: - Status: ok User time: -1.0 System time: 0.04 Wallclock time: 0.05 Termination condition: optimal Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available. Statistics: Branch and bound: Number of bounded subproblems: 0 Number of created subproblems: 0 Black box: Number of iterations: 0 Error rc: 0 Time: 0.12422013282775879 # ---------------------------------------------------------- # Solution Information # ---------------------------------------------------------- Solution: - number of solutions: 0 number of solutions displayed: 0
pyo.value(m.obj)
35.0
m.var_values.extract_values()
{0: 1.0, 1: 1.0, 2: 1.0, 3: 0.0, 4: 1.0, 5: 0.0, 6: 1.0, 7: 1.0, 8: 1.0, 9: 1.0}
m.term_values.extract_values()
{0: 1.0, 1: 0.0, 2: 1.0, 3: 0.0, 4: 1.0, 5: 1.0, 6: 1.0, 7: 1.0, 8: 0.0, 9: 1.0, 10: 0.0, 11: 1.0, 12: 1.0, 13: 1.0, 14: 1.0, 15: 1.0, 16: 0.0, 17: 1.0, 18: 1.0, 19: 1.0, 20: 0.0, 21: 0.0, 22: 0.0, 23: 1.0, 24: 0.0, 25: 1.0, 26: 1.0, 27: 1.0, 28: 1.0, 29: 1.0, 30: 0.0, 31: 1.0, 32: 1.0, 33: 0.0, 34: 1.0, 35: 1.0, 36: 1.0, 37: 0.0, 38: 0.0, 39: 1.0, 40: 1.0, 41: 1.0, 42: 0.0, 43: 1.0, 44: 1.0, 45: 1.0, 46: 0.0, 47: 1.0, 48: 1.0, 49: 1.0}

Вероятностная проверка SAT

raise NotImplementedError("Not yet ready")
--------------------------------------------------------------------------- NotImplementedError Traceback (most recent call last) /var/data/cocalc/1d588dae-0d14-422a-88b4-c470ec2c8303/hard-problems-formulations/minimum-k-satisfiability.ipynb Cell 12 line 1 ----> <a href='vscode-notebook-cell://xn--17-6kce2c.xn--80apqgfe.xn--p1ai/var/data/cocalc/1d588dae-0d14-422a-88b4-c470ec2c8303/hard-problems-formulations/minimum-k-satisfiability.ipynb#X14sdnNjb2RlLXJlbW90ZQ%3D%3D?line=0'>1</a> raise NotImplementedError("Not yet ready") NotImplementedError: Not yet ready
def gen_random_3sat(n_vars: int, n_terms: int) -> np.ndarray['n_terms,3']: """ Returns C - the matrix of the CNF terms. A term is a disjunction of 3 literals. A literal is a var or -var. A var is an integer in [1, n_vars] """ C: np.ndarray['n_terms,3'] = 1 + np.random.randint(n_vars, size=(n_terms, 3)) C *= np.random.choice([-1, 1], size=(n_terms, 3)) return C
gen_random_3sat(6, 10)
array([[ 4, 1, -5], [-1, 5, 4], [-2, 2, 1], [-6, -2, -6], [ 4, -4, -5], [ 2, -5, -1], [ 1, 5, -1], [ 3, 4, -6], [ 4, -4, -1], [ 2, 4, 1]])
def convert_3sat_to_task(C: np.ndarray['n_terms,3']) -> tuple[int, int, int, np.ndarray['n_terms,3']]: k: int = 3 anti_n_terms: int = C.shape[0] n_vars: int = np.max(np.abs(C)) assert C.shape == (n_terms, k) # TODO: ? C = invert_cnf(C) n_terms: int = C.shape[0] assert C.shape == (n_terms, k) return k, n_vars, n_terms, C
def pysat_solve_3sat(C: np.ndarray['n_terms,3']) -> bool: with Solver() as solver: for i in range(C.shape[0]): clause = list(map(int, C[i, ...])) # print(">", clause) solver.add_clause(clause) return solver.solve()
def my_solve_3sat(C: np.ndarray['n_terms,3']) -> bool: m: pyo.ConcreteModel = get_model(*convert_3sat_to_task(C)) solver = pyo.SolverFactory('cbc') solver.solve(m) return np.isclose(pyo.value(m.obj), 0)
def test_3sat(n: int = 100) -> bool: for _ in range(n): C = gen_random_3sat(5, 3) matches: bool = pysat_solve_3sat(C) == my_solve_3sat(C) if not matches: print(f"Discrpeancy found! C:\n{C}") return False return True assert test_3sat() print("Yay!")
Discrpeancy found! C: [[ 2 1 -4] [ 4 3 4] [ 2 2 -4]]
--------------------------------------------------------------------------- AssertionError Traceback (most recent call last) Input In [23], in <cell line: 13>() 9 return False 11 return True ---> 13 assert test_3sat() 14 print("Yay!") AssertionError: